Habitat Selection Tool with Walkthrough for American Pika

Tool Highlights¶

  • This tool can be utilized to determine habitat area based on aspect, slope, and elevation.
  • A walkthrough and appropriate DEM data (downloaded from the USGS) have been provided for identifying American Pika habitat.
  • Alternative DEM data may be utilized and the tool can be customized to desired values of the three included criteria (slope, elevation and aspect).
  • This walkthrough provides instructions on how to customize the tool to your application.
  • Final results can be optionally uploaded to your ArcGIS Online account and visualized on a map widget in the notebook.

Begin Walkthrough¶

Import Modules: The following block of code imports all the modules necessary for the completion of the analysis. The multiprocessingWorker module and corresponding worker function are user defined and where the raster analysis occurs. Its contents (located in the multiprocessing folder of this tool) and how it functions can be viewed by opening it in an IDE or text editor of your choice.

In [ ]:
import time 
import os
import arcpy
import multiprocessingZ
from multiprocessingWorker_v2 import worker

Define Variables: The following variables must be defined or changed by the user for the tool to function properly or to customize the tool to your desired parameters. Comments in the code will help guide you on how to do this successfully.

  • Esri RemapRange Documentation
In [ ]:
'''Point the following variables to the proper location of the test data provided or your own dataset. The sourceWorkspace
is mandatory; however, the clipBoundaryLayer is optional and can be commented out if no boundary is desired for your dataset.'''

sourceWorkspace = r"C:\Users\kevin\Documents\Penn State GIS\Geog 489\Final Project\Test Data_Inyo NF" #Path of the folder that contains all the raster data that you want to process
clipBoundaryLayer = sourceWorkspace + r"\Inyo Boundary\InyoNF_Boundary.shp" #Path to boundary layer for extracting raster to desired area *comment out if no boundary desired

'''The elevCutoff variable requires 3 values in a tuple or list. The first value is a cutoff elevation.
The second is the new value for elevations below the cutoff and third for elevations above the cutoff.'''

elevCutoff = (3500,0,1)

'''The variables below must contain a Remap Table for the desired range of values and a new corresponding
value for the slope and aspect criteria. The predefined tables are for performing the Pika habitat
selection with the provided data.'''

aspectReclassTable = [[-1,-1,0],[0,67.5,1],[292.5,360,1],[67.5,292.5,0]]
slopeReclassTable = [[0,30,0],[30,40,1],[40,90,0]]

Criteria Explanation for Pika Habitat:

  • Elevation: The elevation cutoff splits the raster's elevations into two blocks, those above the cutoff and those below. In this example the cutoff has been set to 3500 meter (~11,500 ft) as this is the appoximate elevation of treeline in our study area (California's Inyo National Forest). Elevations above the cutoff receive a value of "1" as suitable for American Pika and elevations below receive a "0" indicating the area is unsuitable for Pika.
  • Aspect: The aspect table is structed as to assign a value of "1" to northerly aspects (those between 292.5 and 67.5 degrees moving in a clockwise direction) and a "0" value to all other aspects as Pika are most often found in cooler northerly aspects.
  • Slope: Pika are most often found in talus slopes which can form in varying degrees of slope. However, the most common range in which these features will form is between 30-40 degrees of slope. Therefore, we have chosen to assign a value of "1" as suitable for slopes of this degree range in our study are and "0" to all other degree slope ranges.

Set workspace, and populate raster list.

In [ ]:
#Sets arcpy workspace to folder defined in variables block above
arcpy.env.workspace = sourceWorkspace

#Populates a list with all of the rasters that will be sent to workers for geoprocessing
rasterList = arcpy.ListRasters()

Define Multiprocessing Handler

In [ ]:
def mp_handler():
 
    try: 
        
        successList = []
        failedList = []
        jobs = []
        
        '''If no boundary layer used, remove 'clipBoundaryLayer' from jobs.append list'''
        for raster in rasterList:
            jobs.append((raster, sourceWorkspace, elevCutoff, aspectReclassTable, slopeReclassTable, clipBoundaryLayer)) 
    
        cpuNum = multiprocessing.cpu_count()  # determine number of cores to use
        print("There are " + str(cpuNum) + " cpu cores on this machine.")
        print('---------------')
  
        with multiprocessing.Pool(processes=cpuNum) as pool: # Create the pool object 
            returns = pool.starmap(worker, jobs)  # run jobs in job list; returns is a list with return values of the worker function
        
        # Append rasters that were successfully processed and those that failed to process to appropriate list. 
        for res in returns:
            if res.get('result') != False:
                successList.append(res.get('Raster'))
            
            else:
                print(res.get('Raster'))
                print(res.get('error'))
                failedList.append(res.get('Raster'))
         
        #Print success and failed Lists
        if bool(successList) != False:
            print("Successfully processed: " + ', '.join(successList))
        else:
            pass
        
        if bool(failedList) != False:
            print('Failed to process: ' + ', '.join(failedList))
        else:
            print('No rasters failed to process.')
         
    except Exception as e: 
            # Capture all other errors 
            print("Exception:", e)

Execute Analysis Process: Running the block of code below will call on the mp_handler that will generate a list of raster processing jobs that will then be passed to an individual worker. Once complete you can find the final shapefiles in the sourceWorkspace folder designated above. It will also log how long the process takes to complete.

In [ ]:
start_time = time.time() 

if __name__ == '__main__':   
    mp_handler()
    
print ("--- %s seconds ---" % (time.time() - start_time))

Upload Shapefiles to AGOL and Visualize (Optional)¶

Sign in to ArcGIS Online Account

The default requires a Pennsylvania State Organizational account but can be changed to your organization's information

In [ ]:
import arcgis
from arcgis.gis import GIS

#Define your AGOL Username
usernameAGOL = 'kap6389_pennstate' 

#Login to AGOL and paste success code in box that appears below.
gis = GIS('https://pennstate.maps.arcgis.com', client_id='lDSJ3yfux2gkFBYc')
In [ ]:
import zipfile
import re

# function to create zipped shapefile for a given filename without extension
def zipShapefile(name):
    compiledRE = re.compile(name+'(?!.zip)\....')
    with zipfile.ZipFile( os.path.join(sourceWorkspace, name + '.zip'), 'w', zipfile.ZIP_DEFLATED) as zf:  # create zipfile
        for file in os.listdir(sourceWorkspace):                                                           # go through files in workspace
            if compiledRE.match(file):                                                               # test whether file is part of the shapefile to be zipped
                zf.write(os.path.join(sourceWorkspace,file),file,zipfile.ZIP_DEFLATED)                     # add file to zipfile
In [ ]:
# search for content in your gis with a query built from tilte, owner and item type
def searchAGOL(title, owner, itemType):
    return gis.content.search(query='title:'+title+' owner:'+owner, item_type=itemType)

# test whether items exist on AGOL for given title, owner, and item type and if so, delete them from AGOL
def deleteIfExistsOnAGOL(title, owner, itemType):
    result = searchAGOL(title, owner, itemType)   # search item
    print('Found items for title='+title+', owner='+owner+', itemType='+itemType+':')
    print(result)
    for item in result:                           # delete items found
        item.delete()
        print('Item ' + item.title + ' has been deleted.')
In [ ]:
finalMap = gis.map()
finalMap
In [ ]:
shapefileList = arcpy.ListFeatureClasses()
for shapefile in shapefileList:
    name = (os.path.splitext(shapefile)[0])
    zipShapefile(name)
    deleteIfExistsOnAGOL(name, usernameAGOL, 'Shapefile')
    deleteIfExistsOnAGOL(name, usernameAGOL, 'Feature Service')
    Shapefile = gis.content.add({'type':'Shapefile'},  os.path.join(sourceWorkspace, name +'.zip'))     #add shapefile to AGOL
    FeatureService = Shapefile.publish()                                #Publish shapefile as Feature Service
    print(searchAGOL(name, usernameAGOL, 'Shapefile'))                  #Check if Shapefile uploaded
    print(searchAGOL(name, usernameAGOL, 'Feature Service'))            #Check if feature service published
    finalMap.add_layer(FeatureService,{
               'renderer': 'ClassedColorRenderer', 
               'field_name': 'gridcode', 
              })
In [ ]: